Compressive X-ray phase tomography based on transport of intensity 
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We develop a compressive reconstruction method for tomographic recovery of refractive index distribution for 
weakly attenuating objects by using the transport of intensity imaging in a microfocus X-ray system. 



Traditional tomography with hard X-rays recovers the 
attenuation of an object. Attenuation does not always 
provide good contrast when imaging objects made of 
materials with low electron density, e.g. soft tissues. In 
these cases, richer information is often contained in the 
phase, i.e. the optical thickness of the sample [IH3], since 
the X-ray phase shift is almost 10^ times larger than 
the X-ray attenuation for low-Z elements [3]. Propaga- 
tion based techniques are particularly suitable for X-ray 
phase imaging because they allow phase to be recovered 
from intensity images taken at multiple propagation dis- 
tances without the need for optical elements [5l|6]. Here 
we adopt the transport of intensity equation (TIE) which 
relates the measured intensity to the Laplacian of the 
phase under a weakly-attenuating sample approxima- 
tion. Implementing TIE at many angles while rotating 
the object allows tomographic reconstruction of the re- 
fractive index distribution [7l[8]. 

TIE tomographic reconstruction consists of two tasks: 
(1) retrieving phase projections by solving the TIE at 
each angle and (2) applying a tomographic inversion 
algorithm to the phase projections to reconstruct the 
refractive index distribution. The two steps are often 
treated as separate problems |9j. First, the TIE is solved 
by a Poisson equation solver. This step typically re- 
quires regularization, e.g. Tikhonov [9] to mitigate am- 
plification of low spatial frequency noise by the inverse 
Laplacian operator. For the second step, standard tomo- 
graphic reconstruction is carried out, e.g. using the fil- 
tered back-projection (FBP) method. This requires ad- 
ditional regularization to mitigate the well known streak- 
ing (high-frequency) artifacts due to under-sampling in 
the Fourier domain. In this Letter, we design a forward 
model that combines TIE and tomography operations in 
a single, discretized linear operator and develop a com- 
pressive reconstruction method which allows suppression 
of both low and high-frequency artifacts. 

A schematic diagram of the imaging geometry is shown 
in Fig. [H A point source with mean wavelength A and 
intensity /q is located at the plane of z = —zq. This is 
a good approximation for table-top microfocus X-ray 
source. The (weakly attenuating) sample is character- 
ized by a real-valued refractive index n(x,^,z), where 



the origin of the cartesian coordinates is located within 
the object. We further assume that the the object is small 
enough that the beam passing through it can be approx- 
imated as a plane wave oriented along the z axis and 
that the interaction between the sample and the field 
can be treated using the projection approximation, i.e. 
phase delay imparted upon the field passing through the 
sample is ^ = (27r/A) J n(x, y^ z)dz. The intensity / is 
recorded by an area detector located at {x' ^ y' ^d). Under 
the paraxial and small-wavelength approximations, the 
relationship between / and (j) is given by [10] 
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where Vxy is the gradient operator in the {x^y) plane, 
I\ is the background image taken without the sample in 
place, Ms = (zo + d)/ zq is a geometrical magnification 
factor, and d' = z^d/^ZQ + d) is the effective propaga- 
tion distance. Equation [T] is a form of the TIE that uses 
two images taken with and without the sample in place 
in order to recover the projected phase of a weakly at- 
tenuating sample. In a tomographic measurement, the 
sample is rotated about the x axis and the projected 
phase (j) for a given rotation angle can be expressed as 

(l){x,y;0) = n{x,ys,Zs)S{y-ysCose-ZsSmO)dysdZs, 

(2) 
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Fig. 1: Imaging geometry for TIE tomography 



where parametrizes the angle of rotation. The set of 
measurements obtained from many rotation angles can 
be arranged into a vector g. The forward model that 
relates the unknown n to the data takes the simple form 



g = PR n = An. 



(3) 



where the TIE operator P is performed on each phase 
projection and the forward operator A = PR takes the 
projection of n by i? at each angle. 

Since n generally contains more unknowns than 
measurements in g, inversion is ill-posed; assuming n 
can be expressed as sparse, i.e. it contains only a small 
number of nonzero coefficients in some specified basis, 
an effective and elegant approach to invert it is to adapt 
compressed sensing theory [11] . Since the sample of in- 
terest often consists of regions with constant refractive 
indices, we choose the following compressive reconstruc- 
tion model 
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arg min ||n||TV such that g = An, 
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where the total variation (TV) function 
||n||TV is our sparsity basis, defined as 
Ihllxv = E V(V.n)2 + (v^^)2 + (v,n)2, and V„ 
Vy, and \/z are the finite difference operators in the 
three spatial dimensions [12]. We adapt the two-step 
iterative shrinkage/thresholding algorithm (TwIST) [13] 
to solve the minimization. 

When implementing the algorithm, n is discretized 
into di N X N X N matrix with voxel length A. The 
operator P is performed in the Fourier domain: 



^ly = ^xyH oFxy, 



(5) 



where T^y and T~y denote the discrete Fourier and in- 
verse Fourier transform, respectively, with respect to the 
(x, y) variables. H is the transfer function matrix with 
entries defined as Hpq = -47r^(p^ + q'^)/{NA)'^; and o 
denotes entry-wise multiplication. The operator R can 
be implemented in either the spatial domain [14] or the 
Fourier domain [15]. Generally, it is observed that the 
performance of the Fourier domain method is more ro- 
bust to discretization and noise, which has been demon- 
strated in both X-ray CT ^ and MRI [H]. We adopt 
the Fourier domain method to write R as 



projection operators may be combined and simplified as 

^ = JP-^/JoJ-,J-NU, (7) 

where J-nu is the NUFFT operator. 

A microfocus source (Hamamatsu, L8121-03) with a 
circular focal spot size of 5/im in diameter, located at 
zo = 0.765, was operated at 20k Vp to produce a diverg- 
ing X-ray beam with central wavelength A = 0.062nm. 
For a beetle sample, intensity images were taken with a 
custom designed high sensitivity EMC CD based X-ray 
camera, where a 150/im thick RMD CsLTl scintillator 
was fiber optically coupled to an EMCCD (Andor, ixon, 
512x512 pixels, 16/im pixel size) with a 6:1 taper [TS] : 
the effective pixel size at the scintillator is 96/im. The 
scintillator was placed at d = 1.711m. We assume the 
beetle sample is weakly attenuating. During the tomo- 
graphic measurement, a single image was taken at every 
5 degrees with 6.7 seconds exposure time. The intensity 
of the incident beam !{ was calibrated by taking a single 
background image without the sample in place. 

During the reconstruction, we first compute g for each 
angle, an example of which is shown in Fig. [2l Recon- 
struction results using two different methods are com- 
pared in Fig. [3l In (a), the Fourier domain TIE solver 
with Tikhonov regularization chosen to provide optimal 
results is used to compute phase projections at each an- 
gle, and then the FBP method with a Ram-Lak filter is 
applied for the tomographic inversion. The results suffer 
from severe streaking artifacts due to missing samples 
between slices in the Fourier domain. Low-frequency ar- 
tifacts (blurring) around edges are also observable but 
are less severe as compared to a single projected phase 
reconstruction. This is likely due to denser sampling 
around the origin resulting from the intersection of the 
Fourier slices. Both artifacts can be greatly suppressed 
using compressive reconstruction with TV minimization, 
Eq. (J4]), whose results are shown in (b), since TV min- 
imization favors large structures with sharp edges. 3D 
renderings of the compressive reconstruction of the re- 
fractive index are shown in Fig. ^c) and (d). 
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where Tyz,x denotes the discrete 2D Fourier trans- 
form taken for the cross section at a constant x, 
S takes radial slices in the Fourier domain, and 

J^y^ = l^yX^^yX^ ' ' ' ^ ^y}j^ ^^^^^ ^y,L ^enotcs the 
ID inverse Fourier transform of a radial slice along the 
angle Om- Since the radial slices need to be taken from 
points equally spaced along the radial direction, while 
the result of J^yz,x is on a Cartesian grid, S requires 
interpolation and then resampling (gridding) [15]. We 
implement S together with Tyz^x by adapting the non- 
uniform FFT (NUFFT) algorithm [15]. The TIE and 



Fig. 2: Single projection of a beetle sample computed 
according to ([T]). 




Fig. 3: Reconstruction results for the refractive index, (a) Fourier based TIE solver + FBP; (b) Compressive re- 
construction. The three cross-sections are taken from the three orthogonal planes going through the center of the 
sample; (c-d) Two perspectives of a 3D rendering of the reconstruction in (b). [Movies show the rotation of the 
renderings in (c) (Media 1) and (d) (Media 2) about their respective vertical axes.] 



Compressed sensing works best when the measure- 
ment is "incoherent," i.e. the sparse information in the 
unknown is evenly spread out in the measurement [TTi. 
This is achieved in our model by the projection opera- 
tor R. The measurement could be made more incoher- 
ent through the use of source coding [19| or coded aper- 
tures [20] . This is beyond the scope of our current work. 
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